{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 工作年限与收入之间的散点图\n",
    "# 导入第三方模块\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "# 导入数据集\n",
    "income = pd.read_csv(r'C:\\Users\\Administrator\\Desktop\\Salary_Data.csv')\n",
    "# 绘制散点图\n",
    "sns.lmplot(x = 'YearsExperience', y = 'Salary', data = income, ci = None)\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 简单线性回归模型的参数求解\n",
    "# 样本量\n",
    "n = income.shape[0]\n",
    "# 计算自变量、因变量、自变量平方、自变量与因变量乘积的和\n",
    "sum_x = income.YearsExperience.sum()\n",
    "sum_y = income.Salary.sum()\n",
    "sum_x2 = income.YearsExperience.pow(2).sum()\n",
    "xy = income.YearsExperience * income.Salary\n",
    "sum_xy = xy.sum()\n",
    "# 根据公式计算回归模型的参数\n",
    "b = (sum_xy-sum_x*sum_y/n)/(sum_x2-sum_x**2/n)\n",
    "a = income.Salary.mean()-b*income.YearsExperience.mean()\n",
    "# 打印出计算结果\n",
    "print('回归参数a的值：',a)\n",
    "print('回归参数b的值：',b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept          11.346940\n",
       "YearsExperience    24.950094\n",
       "dtype: float64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 导入第三方模块\n",
    "import statsmodels.api as sm\n",
    "# 利用收入数据集，构建回归模型\n",
    "fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit()\n",
    "# 返回模型的参数值\n",
    "fit.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 多元线性回归模型的构建和预测\n",
    "# 导入模块\n",
    "from sklearn import model_selection\n",
    "# 导入数据\n",
    "Profit = pd.read_excel(r'C:\\Users\\Administrator\\Desktop\\Predict to Profit.xlsx')\n",
    "# 将数据集拆分为训练集和测试集\n",
    "train, test = model_selection.train_test_split(Profit, test_size = 0.2, random_state=1234)\n",
    "# 根据train数据集建模\n",
    "model = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + C(State)', data = train).fit()\n",
    "print('模型的偏回归系数分别为：\\n', model.params)\n",
    "# 删除test数据集中的Profit变量，用剩下的自变量进行预测\n",
    "test_X = test.drop(labels = 'Profit', axis = 1)\n",
    "pred = model.predict(exog = test_X)\n",
    "print('对比预测值和实际值的差异：\\n',pd.DataFrame({'Prediction':pred,'Real':test.Profit}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 生成由State变量衍生的哑变量\n",
    "dummies = pd.get_dummies(Profit.State)\n",
    "# 将哑变量与原始数据集水平合并\n",
    "Profit_New = pd.concat([Profit,dummies], axis = 1)\n",
    "# 删除State变量和California变量（因为State变量已被分解为哑变量，New York变量需要作为参照组）\n",
    "Profit_New.drop(labels = ['State','New York'], axis = 1, inplace = True)\n",
    "\n",
    "# 拆分数据集Profit_New\n",
    "train, test = model_selection.train_test_split(Profit_New, test_size = 0.2, random_state=1234)\n",
    "# 建模\n",
    "model2 = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + Florida + California', data = train).fit()\n",
    "print('模型的偏回归系数分别为：\\n', model2.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 导入第三方模块\n",
    "import numpy as np\n",
    "# 计算建模数据中，因变量的均值\n",
    "ybar = train.Profit.mean()\n",
    "# 统计变量个数和观测个数\n",
    "p = model2.df_model\n",
    "n = train.shape[0]\n",
    "# 计算回归离差平方和\n",
    "RSS = np.sum((model2.fittedvalues-ybar) ** 2)\n",
    "# 计算误差平方和\n",
    "ESS = np.sum(model2.resid ** 2)\n",
    "# 计算F统计量的值\n",
    "F = (RSS/p)/(ESS/(n-p-1))\n",
    "print('F统计量的值：',F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model2.fvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 导入模块\n",
    "from scipy.stats import f\n",
    "# 计算F分布的理论值\n",
    "F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1)\n",
    "print('F分布的理论值为：',F_Theroy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 模型的概览信息\n",
    "model2.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 正态性检验\n",
    "# 直方图法\n",
    "# 导入第三方模块\n",
    "import scipy.stats as stats\n",
    "# 中文和负号的正常显示\n",
    "plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "# 绘制直方图\n",
    "sns.distplot(a = Profit_New.Profit, bins = 10, fit = stats.norm, norm_hist = True,\n",
    "             hist_kws = {'color':'steelblue', 'edgecolor':'black'}, \n",
    "             kde_kws = {'color':'black', 'linestyle':'--', 'label':'核密度曲线'}, \n",
    "             fit_kws = {'color':'red', 'linestyle':':', 'label':'正态密度曲线'})\n",
    "# 显示图例\n",
    "plt.legend()\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 残差的正态性检验（PP图和QQ图法）\n",
    "pp_qq_plot = sm.ProbPlot(Profit_New.Profit)\n",
    "# 绘制PP图\n",
    "pp_qq_plot.ppplot(line = '45')\n",
    "plt.title('P-P图')\n",
    "# 绘制QQ图\n",
    "pp_qq_plot.qqplot(line = 'q')\n",
    "plt.title('Q-Q图')\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 导入模块\n",
    "import scipy.stats as stats\n",
    "stats.shapiro(Profit_New.Profit)\n",
    "\n",
    "# 生成正态分布和均匀分布随机数\n",
    "rnorm = np.random.normal(loc = 5, scale=2, size = 10000)\n",
    "runif = np.random.uniform(low = 1, high = 100, size = 10000)\n",
    "# 正态性检验\n",
    "KS_Test1 = stats.kstest(rvs = rnorm, args = (rnorm.mean(), rnorm.std()), cdf = 'norm')\n",
    "KS_Test2 = stats.kstest(rvs = runif, args = (runif.mean(), runif.std()), cdf = 'norm')\n",
    "print(KS_Test1)\n",
    "print(KS_Test2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 导入statsmodels模块中的函数\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "# 自变量X(包含RD_Spend、Marketing_Spend和常数列1)\n",
    "X = sm.add_constant(Profit_New.ix[:,['RD_Spend','Marketing_Spend']])\n",
    "\n",
    "# 构造空的数据框，用于存储VIF值\n",
    "vif = pd.DataFrame()\n",
    "vif[\"features\"] = X.columns\n",
    "vif[\"VIF Factor\"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "# 返回VIF值\n",
    "vif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 计算数据集Profit_New中每个自变量与因变量利润之间的相关系数\n",
    "Profit_New.drop('Profit', axis = 1).corrwith(Profit_New.Profit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 散点图矩阵\n",
    "# 导入模块\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "# 绘制散点图矩阵\n",
    "seaborn.pairplot(Profit_New.ix[:,['RD_Spend','Administration','Marketing_Spend','Profit']])\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 模型修正\n",
    "model3 = smf.ols('Profit ~ RD_Spend + Marketing_Spend', data = train).fit()\n",
    "# 模型回归系数的估计值\n",
    "model3.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 异常值检验\n",
    "outliers = model3.get_influence()\n",
    "\n",
    "# 高杠杆值点（帽子矩阵）\n",
    "leverage = outliers.hat_matrix_diag\n",
    "# dffits值\n",
    "dffits = outliers.dffits[0]\n",
    "# 学生化残差\n",
    "resid_stu = outliers.resid_studentized_external\n",
    "# cook距离\n",
    "cook = outliers.cooks_distance[0]\n",
    "\n",
    "# 合并各种异常值检验的统计量值\n",
    "contat1 = pd.concat([pd.Series(leverage, name = 'leverage'),pd.Series(dffits, name = 'dffits'),\n",
    "                     pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook, name = 'cook')],axis = 1)\n",
    "# 重设train数据的行索引\n",
    "train.index = range(train.shape[0])\n",
    "# 将上面的统计量与train数据集合并\n",
    "profit_outliers = pd.concat([train,contat1], axis = 1)\n",
    "profit_outliers.head()\n",
    "\n",
    "# 计算异常值数量的比例\n",
    "outliers_ratio = sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/ccpp_outliers.shape[0]\n",
    "outliers_ratio\n",
    "# 挑选出非异常的观测点\n",
    "none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu)<=2,]\n",
    "\n",
    "# 应用无异常值的数据集重新建模\n",
    "model4 = smf.ols('Profit ~ RD_Spend + Marketing_Spend', data = none_outliers).fit()\n",
    "model4.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Durbin-Watson统计量\n",
    "# 模型概览\n",
    "model4.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 方差齐性检验\n",
    "# 设置第一张子图的位置\n",
    "ax1 = plt.subplot2grid(shape = (2,1), loc = (0,0))\n",
    "# 绘制散点图\n",
    "ax1.scatter(none_outliers.RD_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())\n",
    "# 添加水平参考线\n",
    "ax1.hlines(y = 0 ,xmin = none_outliers.RD_Spend.min(),xmax = none_outliers.RD_Spend.max(), color = 'red', linestyles = '--')\n",
    "# 添加x轴和y轴标签\n",
    "ax1.set_xlabel('RD_Spend')\n",
    "ax1.set_ylabel('Std_Residual')\n",
    "\n",
    "# 设置第二张子图的位置\n",
    "ax2 = plt.subplot2grid(shape = (2,1), loc = (1,0))\n",
    "# 绘制散点图\n",
    "ax2.scatter(none_outliers.Marketing_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())\n",
    "# 添加水平参考线\n",
    "ax2.hlines(y = 0 ,xmin = none_outliers.Marketing_Spend.min(),xmax = none_outliers.Marketing_Spend.max(), color = 'red', linestyles = '--')\n",
    "# 添加x轴和y轴标签\n",
    "ax2.set_xlabel('Marketing_Spend')\n",
    "ax2.set_ylabel('Std_Residual')\n",
    "\n",
    "# 调整子图之间的水平间距和高度间距\n",
    "plt.subplots_adjust(hspace=0.6, wspace=0.3)\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# BP检验\n",
    "sm.stats.diagnostic.het_breushpagan(model4.resid, exog_het = model4.model.exog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 模型预测\n",
    "# model4对测试集的预测\n",
    "pred4 = model4.predict(exog = test.ix[:,['RD_Spend','Marketing_Spend']])\n",
    "# 对比预测值与实际值\n",
    "pd.DataFrame({'Prediction':pred4,'Real':test.Profit})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 绘制预测值与实际值的散点图\n",
    "plt.scatter(x = test.Profit, y = pred4)\n",
    "# 添加斜率为1，截距项为0的参考线\n",
    "plt.plot([test.Profit.min(),test.Profit.max()],[test.Profit.min(),test.Profit.max()],\n",
    "        color = 'red', linestyle = '--')\n",
    "# 添加轴标签\n",
    "plt.xlabel('实际值')\n",
    "plt.ylabel('预测值')\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
